subroutine solver
    use datatype
    use omp_lib
    implicit none

    write (*, *)
    write (*, '(A19)') 'CFD process begin !'
    do i = 1, n
        do j = 1, m
            allocate (cell(i, j)%phi(1:2))
            cell(i, j)%u = 0.0
            cell(i, j)%v = 0.0
            cell(i, j)%cp = 0.0
            cell(i, j)%phi(:) = 0.0
        end do
    end do

    do i = 1, n
        do j = 1, m
            cell(i, j)%phi(1:2) = v_inf*cell(i, j)%x(2) !< init
        end do
    end do

    open (unit=123, file='cp.plt', form='formatted')
    !$omp parallel num_threads(numprocs)

    do iter = 1, maxstp_cfd !< loop

        !< bc: far field & periodic
        do i = 1, m
            do j = 1, n
                if (cell(j, i)%y(1) == 0.0 .and. cell(i, j)%x(1) >= c) then !< left & right: periodic
                    cell(j, 1)%phi(1) = 0.5*(cell(j, 2)%phi(1) + cell(j, m - 1)%phi(1))
                    jacobian(j, 1) = 0.5*(jacobian(j, 2) + jacobian(j, m - 1))
                    cell(j, m)%phi(1) = cell(j, 1)%phi(1)
                end if
                cell(n, i)%phi(1) = v_inf*cell(n, 1)%x(2)!< top: farfield
            end do
        end do

        !< bc: solid wall
        do j = 2, m - 1 !< special work for inner-wall cells
            cell(1, j)%phi(2) = cell(2, j)%phi(1) - 0.5*beta(1, j)/gamma(1, j)* &
                                (cell(1, j + 1)%phi(1) - cell(1, j - 1)%phi(1))

            jacobian(1, j) = (cell(1, j + 1)%x(2) - cell(1, j - 1)%x(2))*(cell(2, j)%y(2) - cell(1, j)%y(2))/2 - &
                             (cell(1, j + 1)%y(2) - cell(1, j - 1)%y(2))*(cell(2, j)%x(2) - cell(1, j)%x(2))/2
            jacobian(n, j) = (cell(n, j + 1)%x(2) - cell(n, j - 1)%x(2))*(cell(n - 1, j)%y(2) - cell(n, j)%y(2))/2 - &
                             (cell(n, j + 1)%y(2) - cell(n, j - 1)%y(2))*(cell(n - 1, j)%x(2) - cell(n, j)%x(2))/2

            jacobian(1, j) = 1/jacobian(1, j)
            jacobian(n, j) = 1/jacobian(n, j)

            cell(1, j)%u = jacobian(1, j)*((cell(1, j + 1)%phi(1) - cell(1, j - 1)%phi(1))* &
                                           (cell(2, j)%y(2) - cell(1, j)%y(2)) - &
                                           (cell(2, j)%phi(1) - cell(1, j)%phi(1))*(cell(1, j + 1)%y(2) - cell(1, j - 1)%y(2)))/2

            cell(1, j)%v = jacobian(1, j)*((cell(2, j)%phi(1) - cell(1, j)%phi(1))* &
                                           (cell(1, j + 1)%x(2) - cell(1, j - 1)%x(2)) - &
                                           (cell(1, j + 1)%phi(1) - cell(1, j - 1)%phi(1))*(cell(2, j)%x(2) - cell(1, j)%x(2)))/2

            cell(1, j)%cp = 1 - (cell(1, j)%u**2 + cell(1, j)%v**2)/v_inf**2
        end do

        !< -----------------------------------------------------------------------------------------------------
        cell(1, 1)%phi(2) = 1/(1 - (beta(1, 1)/gamma(1, 1)))*(cell(2, 1)%phi(1) - (beta(1, 1)/gamma(1, 1))* &
                                                              cell(1, 2)%phi(1))
        jacobian(1, 1) = (cell(1, 2)%x(2) - cell(1, 1)%x(2))*(cell(2, 1)%y(2) - cell(1, 1)%y(2)) - &
                         (cell(1, 2)%y(2) - cell(1, 1)%y(1))*(cell(2, 1)%x(2) - cell(1, 1)%x(2))
        jacobian(1, 1) = 1/jacobian(1, 1)
        !< ------------------------------------------------------------------------------------------------------
        cell(1, m)%phi(2) = 1/(1 - (beta(1, m)/gamma(1, m)))*(cell(2, m)%phi(1) + (beta(1, m)/gamma(1, m))* &
                                                              cell(1, m - 1)%phi(1))
        jacobian(1, m) = (cell(1, m)%x(2) - cell(1, m - 1)%x(2))*(cell(2, m)%y(2) - cell(1, m)%y(2)) - &
                         (cell(1, m)%y(2) - cell(1, m - 1)%y(1))*(cell(2, m)%x(2) - cell(1, m)%x(2))
        jacobian(1, m) = 1/jacobian(1, m)
        !< ------------------------------------------------------------------------------------------------------
        jacobian(n, 1) = (cell(n, 2)%x(2) - cell(n, 1)%x(2))*(cell(n - 1, 1)%y(2) - cell(n, 1)%y(2)) - &
                         (cell(n, 2)%y(2) - cell(n, 1)%y(1))*(cell(n - 1, 1)%x(2) - cell(n, 1)%x(2))
        jacobian(n, m) = (cell(n, m)%x(2) - cell(n - 1, m)%x(2))*(cell(n, m)%y(2) - cell(n, m - 1)%y(2))/2 - &
                         (cell(n, m)%y(2) - cell(n, m - 1)%y(2))*(cell(n, m)%x(2) - cell(n - 1, m)%x(2))/2
        jacobian(n, 1) = 1/jacobian(n, 1)
        jacobian(n, m) = 1/jacobian(n, m)

        !< inner cells
        !$omp do
        do i = 2, n - 1
            do j = 2, m - 1

                cell(i, j)%phi(2) = 0.5*(alpha(i, j)*(cell(i, j + 1)%phi(1) + cell(i, j - 1)%phi(1)) + &
                                         gamma(i, j)*(cell(i + 1, j)%phi(1) + cell(i - 1, j)%phi(1)) - &
                               0.5*beta(i, j)*(cell(i + 1, j + 1)%phi(1) + cell(i - 1, j - 1)%phi(1) - cell(i - 1, j + 1)%phi(1) - &
                                                         cell(i + 1, j - 1)%phi(1)))/(alpha(i, j) + gamma(i, j))

                cell(i, j)%res(3) = abs(cell(i, j)%phi(2) - cell(i, j)%phi(1))
                error(i, j, 3) = cell(i, j)%phi(2) - cell(i, j)%phi(1)

                jacobian(i, j) = (cell(i, j + 1)%x(2) - cell(i, j - 1)%x(2))*(cell(i + 1, j)%y(2) - cell(i - 1, j)%y(2))/4 - &
                                 (cell(i, j + 1)%y(2) - cell(i, j - 1)%y(2))*(cell(i + 1, j)%x(2) - cell(i - 1, j)%x(2))/4

                jacobian(i, j) = 1/jacobian(i, j)

                !cell(i, j)%u = (cell(i, j +1)%phi(2) - cell(i, j - 1)%phi(2))/(cell(i, j +1)%x(2) - cell(i, j - 1)%x(2))
                !cell(i, j)%v = (cell(i +1, j)%phi(2) - cell(i - 1, j)%phi(2))/(cell(i + 1, j)%y(2) - cell(i - 1, j)%y(2))

                cell(i, j)%u = jacobian(i, j)*((cell(i, j + 1)%phi(1) - cell(i, j - 1)%phi(1))* &
                                               (cell(i + 1, j)%y(2) - cell(i - 1, j)%y(2)) - &
                                      (cell(i + 1, j)%phi(1) - cell(i - 1, j)%phi(1))*(cell(i, j + 1)%y(2) - cell(i, j - 1)%y(2)))/4

                cell(i, j)%v = jacobian(i, j)*((cell(i + 1, j)%phi(1) - cell(i - 1, j)%phi(1))* &
                                               (cell(i, j + 1)%x(2) - cell(i, j - 1)%x(2)) - &
                                      (cell(i, j + 1)%phi(1) - cell(i, j - 1)%phi(1))*(cell(i + 1, j)%x(2) - cell(i - 1, j)%x(2)))/4

                cell(i, j)%cp = 1 - (cell(i, j)%u**2 + cell(i, j)%v**2)/v_inf**2

            end do
        end do
        !$omp end do
        do i = 1, n
            do j = 2, m
                cell(i, j)%phi(1) = cell(i, j)%phi(2) !< update
            end do
        end do
        max_res(3) = abs(maxval(error(2:nn, 2:mm, 3)))
!$      end_time = omp_get_wtime()

        !$omp barrier
        !$omp single
        if (mod(iter, 5000) == 0) then
            write (*, '(A40)') 'Max error(phi)     Iters     Time'
        end if
        if (mod(iter, 500) == 0) then
            write (*, '(E20.12, I11, F10.2)') max_res(3), iter - 1, end_time - start_time
        end if

        !$omp end single

    end do
    !$omp end parallel

    write (*, 130)


    write (123, *) 'variables = "x", "-Cp"'
    do j = 1, m
        write (123, *) cell(1, j)%x(2), cell(1, j)%cp
    end do

130 format(//9x, '--------- All computation is finished!---------'//)

    return
end subroutine
